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APPLICATION OF LAX'S FINITE-DIFFERENCE METHOD TO 
NONEQUILIBRIUM HYPERSONIC FLOW PROBLEMS 

By Fred R. DeJarnette 
Langley Research Center 

SUMMARY ^ /O *¥ 

The finite-difference method developed by P. D. Lax for unsteady flows is applied 
to the steady- state equations for the supersonic region of inviscid flows past two- 
dimensional and axisymmetric bodies. A diatomic gas subject to nonequilibrium vibra- 
tions and dissociation is considered. As a consequence of Lax's finite-difference scheme, 
an artificial viscosity is implicitly introduced. This scheme allows the computations to 
proceed downstream of an initial data line as if no shock wave were present at all. The 
shock wave appears in the solution, however, smeared over several mesh spaces while it 
accurately gives the proper jump conditions across the shock. 

A new difference scheme is developed for the body points that is more accurate 
than previously used extrapolation schemes and reflection principles. The numerical 
results of the present method are shown to compare favorably with the methods of char- 
acteristics and integral relations for frozen flow over a cone — parabolic -arc— cylinder, 
vibrational nonequilibrium flow over a wedge, and chemical nonequilibrium flow of 
Lighthill's "ideal dissociating gas" over a wedge. 

INTRODUCTION 

The numerical solution of the inviscid flow field surrounding a vehicle traveling at 
hypersonic speeds provides a formidable task for computation. At the speeds encountered 
for a planetary entry, the problem is made more difficult by the presence of nonequilib- 
rium real gas effects. This report develops a finite-difference method for computing the 
supersonic-hypersonic regions of steady, inviscid flows past two-dimensional and axisym- 
metric bodies. 

Many methods have been formulated and used for nonequilibrium flow-field compu- 
tations varying from the accurate, but lengthy, method of characteristics (refs. 1 to 3) and 
inverse blunt-body methods (ref. 4) to the gross, yet simple, method of integral relations 
(refs. 5 to 7) and stream-tube techniques (ref. 8). In 1954 Lax (ref. 9) introduced a first- 
order-accuracy difference scheme for computing one-dimensional unsteady flow fields 
containing shock waves. Roberts (ref. 10) applied this method to spherical waves, and 




4 


Bohachevsky et al. (ref. 11) used it to compute steady- state flow fields for two and three 
dimensions by obtaining the asymptotic solution of the time-dependent flow. More 
recently, Lax and Wendroff (ref. 12) proposed a finite-difference scheme that is accurate 
to the second order. Burstein (ref. 13) used this technique to obtain steady flow fields as 
the asymptotic form of the unsteady flow, whereas Thommen and D’Attorre (ref. 14) 
applied the difference scheme to the steady-flow equations. 

Although the Lax-Wendroff method is accurate to the second order, the numerical 
computations become prohibitive for a real gas. The numerical method presented here 
applies the basic difference scheme of the earlier Lax method (ref. 9) to the steady-state 
flow-field equations. This method represents a good compromise between accuracy and 
numerical complexities. The governing partial differential equations are first written in 
conservation (or divergence) form; that is, the coefficients of the derivatives are all 
unity. The partial derivatives are then replaced with finite-difference quotients, but in 
such a manner as to introduce implicitly an artificial viscosity. However, the artificial 
viscosity is due only to the finite-difference scheme employed and differs from that used 
by VonNeumann and Richtmyer (ref. 15). (VonNeumann and Richtmyer add a viscous 
pressure term to the static pressure in the differential equations before the difference 
equations are formed.) As a consequence of the artificial viscosity effect, the computa- 
tions may proceed downstream of an initial data line as if no shock wave were present at 
all. The shock wave appears in the solution, however, not as a discontinuity but as a near 
discontinuity smeared over several mesh spaces while the solution accurately gives the 
proper jump conditions across the shock wave. 

In order to obtain good accuracy on the body surface, a different difference scheme 
is developed for the body points that is more accurate than the previously used extrapola- 
tion schemes (ref. 13) and reflection principles (refs. 11 and 13). The numerical results 
of the present method are compared with those of the methods of characteristics and 
integral relations for frozen flow over a cone — parabolic-arc — cylinder, vibrational non- 
equilibrium flow over a wedge, and chemical nonequilibrium flow of Lighthill's "ideal 
dissociating gas" over a wedge. The method presented herein is developed in more detail 
in reference 16. 


SYMBOLS 

a m vector defined by equation (16) 

b m vector defined by equation (16) 

C* constant in dissociational rate equation 
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C F locally frozen speed of sound 

* , 2(p - p^) 

C n pressure coefficient, 22Z 

y p V 2 

D dissociation energy per molecule 

d m vector defined by equation (16) 

E v vibrational energy per unit mass 

F function defined by equation (C2) 

G function defined by equation (C3) 

H static enthalpy per unit mass 

Hj. total enthalpy per unit mass 

j integer in flow equation; 0 for two-dimensional flow, 1 for axially symmetric 

flow 

K body curvature 

kg Boltzmann constant 

L length scale 

M locally frozen Mach number 

n,k number of increments in mesh system (see fig. 2) 

N 0 Avogadro's number 

p pressure 

r perpendicular distance from body axis 

R universal gas constant 

t time 

T translational temperature 

u velocity in x-direction 

V magnitude of total velocity 
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V 

W 

X 

y 

a 

r 


y 

6 


©v 

e 

A 


M 

€ 

P 

P D 

T 

* 

Wj,W2> w 3> Ct ’ 

n 


velocity in y-direction 

molecular weight of jth species 

distance along body surface from leading edge! 

> (see fig. 1) 

distance normal to body surface J 

atom mass fraction 

rate of dissociation 

locally frozen ratio of specific heats 

y-position of shock wave 

characteristic dissociation temperature, ©/kg 

characteristic vibrational temperature 

angle between body tangent and axis of symmetry 

vibrational energy rate 

scale factor in x-direction 

locally frozen Mach angle 

Cartesian coordinate along body axis 

mass density 

characteristic dissociation density 
vibrational relaxation time 
constant defined by equation (10) 

functions defined by equations (41) and (47) 
shock-wave angle 


Subscripts: 

A 2 diatomic species 

e equilibrium conditions 
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F 


locally frozen value 


n,k spatial position, x = x Q + nix, y = kAy (see fig. 2) 

o initial conditions 

oo undisturbed free-stream conditions 

Barred quantities are dimensionless as given by equation (7). Unbarred quantities 
are dimensional , 


PROBLEM DEFINITION 
General Description 

The problem to be studied is the numerical solution of the inviscid, steady super- 
sonic flow over two-dimensional and axisymmetric bodies. The gas model considered is 
a homonuclear, diatomic gas subject to nonequilibrium vibrations and dissociation. 
Translation and molecular rotation modes are considered to be completely excited, and 
electronic excitation and ionization are neglected. 

The numerical technique for the flow-field solution is to replace the governing par- 
tial differential equations with finite-difference equations. Then, with known flow prop- 
erties along an initial data line in the supersonic region, the downstream flow field is 
computed by marching forward in the 
streamwise direction. Since the numer- 
ical technique is restricted to hyperbolic 
equations, the flow field computed must 
remain supersonic. 


Basic Equations 

Figure 1 illustrates the geometry and 
coordinate system. For this body- oriented 
coordinate system, x is measured along 
the surface of the body, y is normal to sur- 
face, and the velocity components in these 
respective directions are u and v. The 
basic flow equations are written in conserva- 
tion (or divergence) form as (see ref. 16): 
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Continuity: 


9(pufj) | a(pvr^x) _ 9 

ax ay 


x- momentum: 


aK p + pu2)f j] d(pvuri\ x ) — . 

; J + -±- ^ + KpuvrJ - jpA v sin 0 = 0 

ax ay x • 


y-momentum: 


afay^^lp.pv^fixj _ + „ 2)fj _ 

3x 3y 


jpA x cos 0*0 


( 1 ) 


( 2 ) 


(3) 


Vibrational energy rate: 

a(pufiE v ) a(pvfjx x E v ) X x rjpLW A2 A 

PH r\~ tt m 0 


ax 


3y 


VooT^R 


Rate of dissociation: 


d(pu^aj a(pvfj A x cnj A x fjpLr ^ 

+ 9 y v„ =0 


3x 


Energy: 


_ -2 -2 _ 

H + u ± v = H t = Constant 


(4) 


(5) 


( 6 ) 


where K = - — is the body curvature, A x = 1 + Ky is the scale factor in the 
dx x 

x-direction, and j = 0 for two-dimensional body and j = 1 for axisymmetric body. 

The nondimensional vibrational energy is E v , a is the atom mass fraction, A is the 

dimensional vibrational rate and F is the dimensional dissociational rate 

dt ’ dt 

All quantities are nondimensionalized as follows (barred quantities are nondimensional): 

H_ 

T 2 


u or v 
u or v = — — 


H = JL 


T = ^ 

Too 


E 


: v w A2 


Ey RT. 


OO OO 


0 D or 0 v = 


p = -£- 

Poo 


_ e D or 0 V 




- x, y, or r 

x, y, or r = i ^ L 


( 7 ) 
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The enthalpy and equation of state are, respectively, 


H = i p 


7 + (1 - o)E y + 

dt 



where 


p = pT(l + a)\p 



The locally frozen ratio of specific heats is 

v - 1 ± 3« 

r 5 + a 


and the locally frozen Mach number is 


M = 



( U 2 + v 2 )wa 2 

y(l + Qf)RT 


1/2 


( 8 ) 

( 9 ) 

( 10 ) 


(ID 


( 12 ) 


Vibrational Nonequilibrium Gas Model 

The vibrational nonequilibrium gas model is for a diatomic gas with no dissociation. 
Therefore, a= 0 everywhere and the dissociational rate equation (eq. (5)) is not needed. 
The vibrational energy rate and relaxation time used here are the same as those used in 
references 1 and 5; the vibrational energy rate is 

a = E v? e - _El 

T 

where the local equilibrium vibrational energy is 


®v,e 

and the vibrational relaxation time r 


RQ V 

WA 2 (e 6v/T - l) 
is given by 


rp oc exp 


1.3965 x 10 6 \ 1//3 


where the temperature is in °K. 
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Lighthill Gas Model 


A convenient simplification is achieved when the gas model is reduced to Lighthill' s 
ideal dissociating gas of reference 17. This gas model is used in this report for compar- 
ative purposes only. The essential features of the gas model are: 


(1) The vibrational energy always has one-half of its fully excited classical value; 
that is, 



(2) The characteristic density of dissociation is considered to be constant. 

As a consequence of these simplifications, equation (4) is not required, the enthalpy 
given by equation (8) is reduced (for Lighthill gas) to 


H = <//[(4 + a)T + a'OjjJ (13) 

and the locally frozen ratio of specific heats given by equation (11) is changed (for 
Lighthill gas) to 

y = 4-p (14) 

With the exception of equations (4), (8), and (11), the flow-field equations given by equa- 
tions (1) to (12) also hold for the Lighthill gas. 

The dissociational rate equation 


T = C*pT' 2 ‘ 5 j(l - o)e" 0 D/ r - ^_o£| 

which was used in references 3 and 7 (where C* is a constant) is also used herein. 

Boupdary Conditions 

Since the velocity at the boundary must be tangent to the body surface, v = 0 at 
y = 0. The numerical technique developed here does not consider the shock wave as a 
boundary or a discontinuity. Therefore, there are no boundary conditions to be satisfied 
at the shock wave since the initial data line extends into the free stream. 

METHOD FOR NUMERICAL COMPUTATIONS 

Let the flow field be divided into an orthogonal mesh system as shown in figure 2. 

If the mesh spacings Ax and Ay remain constant, the coordinates of the mesh point 
n,k are x * x 0 + n Ax and y = k Ay. Note, however, that the arc length from n,k to 
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n+l,k is A x Ax, whereas the linear length from n,k to n,k+l is Ay. (See fig. 3.) 
In some applications Ax may vary with n. For these cases 

n 

x = x 0 + ^T(Ax)i 
i=l 


The numerical technique for the flow-field solution is to replace the differential equations 
with finite-difference equations. Then, with known initial data along the mesh points at 
n = 0 (obtained from some other technique), new data are computed along the mesh points 
at n = 1. This technique is continued for increasing values of n. In this manner the 
flow-field properties are computed at the various mesh points throughout the supersonic 
portion of the flow field. Since the differential equations are hyperbolic, the ratio Ax/Ay 
is bounded. If this bound is exceeded, the numerical computations become unstable. The 
proper stability criterion is discussed subsequently. 

The system of partial differential equations, given by equations (1) to (5), may be 
written as one vector equation: 


9a m 

9x 


+ 


9b m 

9y 


+ d. 


m 


= 0 


(m = 1, 2, 3, 4, 5) (15) 
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where 


3 

l 4 


jour) 

bj = pvrix x 

dj = 0 ^ 

(p + pu2)f J 

b 2 = b x u 

d 2 = Ka g - jpX x sin 6 

a x v 

b 3 = (p + pv2)?3x x 

d 3 = -Ka 2 - jpX x cos 6 

a l E v 

b 4 = b x E v 

-X x fiApLWA2 

d 4 “ V T R 

T oo oo 

a^o! 

b 5 = b l0t 

. -X x r i LpT 

5 y 


(16) 


9a 


m 


9x 


The difference scheme, suggested by Lax (ref. 9), replaces the partial derivative 


n,k 


at the mesh point n,k with the modified forward difference quotient 
( am )n+l,k 2^ am )n,k+l + ( am )n,k-l] 


Ax 


and replaces 


9b 


m 


/n,k 


with the symmetric difference quotient 

( bm )n,k+l - ( b m)n,k-l 
2 Ay 


Since only the mesh points n,k+l and n,k-l are involved at n for both of these dif- 
ference quotients, the term (d m ) n k is replaced by the average value 


2 ( dm )n,k+l ' ( d m) n ,k-l] 

With these substitutions the system of partial differential equations 

(9am\ / 9b r 


*-m 

9x 


J m 


n,k \ 3y / n ,k 


+ ( dm )n,k " 0 


(m = 1, 2, 3, 4, 5) 


are replaced with 


( a m) n +l,k “ 2 ( am )n,k+l + ( a m) n ,k-l “ ( b m) n ,k+l ' ( b m) n ,k-l 



(17) 


for m = 1, 2, 3, 4, 5. In this manner |a m j n+1 k may be evaluated in terms of a m , 
b m , and d m at the mesh points n,k+l and n,k-l. Therefore, if the quantities a m 
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are known at n = 0 and k = 0, 2, 4, . . the values of a m can be determined at all 
points of the staggered mesh n = 1, 2, 3, . . . and k = 1, 3, 5, 7, . . . when n is odd; 
and k = 0, 2, 4, 6 when n is even. (See fig. 2.) 

If equation (17) is used for k = 0, then imaginary flow properties are required at 
the mesh point n,-l which is inside the body. Since the computational procedure cannot 
determine these properties at k = -1, a different finite-difference scheme is developed 
later for the body points. 


The significance of the finite -difference scheme employed here is more readily seen 
when equation (17) is rearranged as 


( a m)n+l,k " ( am )n,k . ( b m) n ,k+l ' ( b m) n ,k-l . u \ 

~T = + TT= + i a mj n 

Ax 2 Ay v /u > 




2 Ax 


( a m ~ ASd m) llt k +1 ~ 2 ( a m ~ Axd m) ntk + ( a m - Axd m) ntk _i 


(Ay) 2 


(18) 


The conventional manner of forming finite differences is to replace ( jP -) with the 

\ 9x / n k 

forward difference quotient ’ 


( a m) n+ i i k ~ ( a m) n ,k 
Ax 


ab 


to replace I — m i with the symmetric difference quotient 

9 y/n,k 

( bm )n,k+l " ( bm )n,k-l 
2 Ay 


and to replace the second derivative 


9 ^( a m ~ Axd m ) 


°y 2 Jn,k 


with 


( a m - ^ d rti) n ,k+l " 2 ( am " ^^in.k + ( a m ~ ASd m) n ,k-l 


(Ay) 2 


Therefore, the differential equation that corresponds to the difference equation (17) using 
conventional difference quotients is 


3a m , 9 t>m , ^ _ (Ay)^ 3^(a m - Axd m ) / 1n \ 

ir + ir m 2 ^s> 9j2 1 ’ 

This differential equation differs from that shown as equation (15) by the addition of the 
second derivative term on the right-hand side. In comparison with the momentum equa- 
tion for a viscous fluid (see ref. 18), the additional term in equation (19) is similar to part 
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of the viscous terms in the viscous momentum equation. However, the "artificial" coef- 
ficient of vis cosity is not a function of the fluid properties, as it is in the viscous 

2 Ax 

momentum equation, but is a function of the mesh spacing only. Also, as the mesh 
spacing approaches zero, this coefficient approaches zero. The effect of the additional 
term involving the second derivative is designated in references 9 and 10 as an "artificial 
viscosity." It tends to smooth out discontinuous and rapidly varying flow variables when 
they arise. 

One such discontinuity that appears in the inviscid solution of supersonic flow over 
bodies is a shock wave. Normally, the shock wave is taken as a boundary, and the 
Rankine-Hugoniot equations specify the boundary conditions. However, a numerical solu- 
tion of the viscous Navier-Stokes equations yields a shock-wave structure in which the 
fluid properties are not discontinuous, but vary continuously over a narrow distance. 

(See ref. 19.) The introduction of "artificial" viscosity in equation (18) gives this same 
qualitative effect, but the shock-wave thickness predicted is too large because the viscos- 
ity effect is not real. A simple example of the numerical computation of the flow of a 
perfect gas through an oblique shock wave illustrates this point and is given in appendix A. 

An important feature of the present method is that shock waves are not regarded as 
discontinuities and the numerical solution proceeds as if no such waves were present at 
all. However, the shock waves appear in. the final solution not as a discontinuity but as 
a near-discontinuity smeared over several mesh spaces. As a result, there are no 
boundary conditions to be satisfied at the shock wave itself; for some fluid problems this 
simplification is a significant one. Although the shock-wave thickness in the solutions 
is incorrect, the flow fields on both sides are accurately described. (See the example 
in appendix A.) The success of the difference scheme presented here is due, in part, to 
the divergence form of the differential equations and, in part, to the artificial viscosity 
effect. The difference scheme of equation (17) is accurate to 0(Ax2). (See ref. 16 for 
the details.) 


Stability Criterion 

Consider the two adjacent mesh points n,k+l and n,k-l in a flow field. As has 
been shown, the flow properties at these two points are used to determine those at the 
point n+l,k. Since the numerical technique used here is restricted to supersonic flow, 
the mesh point n+l,k must lie within the Mach triangle about the points n,k+l and 
n,k-l. (See fig. 3.) This criterion is the classical Courant- Friedrichs- Lewy stability 
criterion (ref. 20) used by Lax (ref. 9) and Roberts (ref. 10). 

The Mach triangle is the triangle formed by the line joining n,k-l and n,k+l, the 
left- running Mach line from n,k-l, and the right- running Mach line from n,k+l. As 
pointed out by Wood and Kirkwood (ref. 21), the appropriate Mach line for a reacting gas 
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Pypa Region permissible for mesh point 
n + l, k by stability criterion 
□ Mesh points already computed 
o New mesh point to be computed 

Figure 3.- Stability criterion. 


is that determined by the locally frozen Mach number; that is, 

sin ^ = — (Mg 1) (20) 

M 

where p is the angle between the locally frozen Mach line and the streamline. There- 
fore, in the x,y coordinate system used here, the slopes of the frozen Mach lines are 
given by 



= tan 


tan' 



± P 


After considerable rearrangement, this equation becomes 


where 



uv 



/ {j2 + y2 

Uf 2 ) 

ii 2 


C F 


2 


yp 

P 


( 21 ) 


13 


m =. 

dx/_ 


Cp 

oo and 


When =J-= 1, f^j = ±oo for v = 0, = +<*> and 


dx ( 

d£ 

dx 


dx 


dy 

vdx/_ 


is finite for v > 0, and 


is finite for v < 0. Therefore, at least one of the Mach lines is 


parallel to the y-axis when = 1. This same result may also be deduced from the 

Cp 

fact that the velocity component normal to a frozen Mach line has the magnitude of the 
frozen speed of sound. 


Equation (21) may be used to write the stability criterion as 


Ax 

Ay 



( 22 ) 


Therefore, in order to have Ax/Ay finite and positive, it is necessary to have 

=H-> 1 

Cp 

Note that since ~~ ^ M, =— >1 is a more stringent requirement than M > 1 . 
Cp Cp 


Individual Fluid Properties 

For mesh points away from the body, k^l, equation (17) gives (a m ) n+1 k in 
terms of the known properties at n,k+l and n,k-l. A method is given below whereby 
the individual fluid properties can be obtained from a m . Since 


it is immediately evident that 


a^ = pur^ 
a 2 = (P + P*i2)fJ 


a 3 = a l v 


a 4 = a l E v 


a 5 = a l 0! 



( 23 ) 
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(24) 


E - a 4 
E v-a7 


a l 


The other properties are more involved. The ratio ^2. gives 

a l pu 
or 

S 2- J2u + E=0 

a l ^ 

The £ term in this equation can be obtained as a function of u, 
P 

the energy equation and the equation of state as 


v2 


P 


H t - (1 - a)E - a0 D ^ - ^ 


u2 


7 + 3a 


2(1 + a) 

Substitute equation (28) into equation (27) to get 


7 + 3q 
1 + a 


v2 


6 + 2a-? a 2 - H t - » - “I®'* - o®D^ - V 

a ' - u* - — u + 


(25) 


v, E 


v> 


(26) 

(27) 

and a from 

(28) 


7 + 3a 


7 + 3a 
2(1 + a) 


= 0 


(29) 


and, the solution of this quadratic equation is 


u = 



8(1 + a)(6 + 2a) 
(7 + 3a) 2 


- (1 - a)E v ^/ - q0 D ^/ - 


2(6 + 2a) 
7 + 3a 



(30) 


From a physical standpoint u must be real and single valued at every position in the 
flow field. Thus, the term under the radical must be positive or zero, and the sign pre- 
ceding this term must be appropriately chosen. To investigate these requirements, the 
restrictions necessary to make equation (30) an identity will be determined. 

From equation (26), 


a 2 p - 

s-= + u 

d l pu 
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and equation (28) yields 


8(1 + o;)(6 + 2a) 


(7 + 3a) 2 


H t - (1 - a) E v + - 4] - 


Hence, the term under the radical in equation (30) becomes 


a 2^ 2 
*1 


8(1 + 

(7 + 3a) 


- (1 . „)!> - ^ - f] . [M] 


x _ (7 + 3a) p 
(5 + a)pu 2 


This relation shows that the term under the square root sign is always nonnegative. 
Since the nondimensional frozen speed of sound is 


Cp 



'zl>\ 1//2 = (l + 3 a 
\p J \5 + a 



equation (30) becomes 


u 


u = 


2 


5 + a \ C F 
(l + 3a/ u2 


+ 1 



4(3 + a) 
7 + 3a 


If it is noted that 


u > 0 

Oia^l 


and 


- vfl 

1/2 

1- Cf2 

52 /J 


u 2 


equation (30) becomes an identity if the positive sign is chosen when 

, C F 2 


0 


and the negative sign is selected for 


u 


C 2 


u 4 


^ 0 
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> 0. There- 


The stability criterion for supersonic flow (eq. (22)) requires that 1 - 
fore, the positive sign must be used in equation (30) and, as a result, 



u 



8(1 + a )(6 + 2 a ) 
(7 + 3a) 2 



2(6 + 2a) 
7 + 3a 


(31) 


The remaining flow properties may now be computed from equations (16) and (9) as 


P = 


ufi 


p = 


a 2 " a l u 


and 


T = 


p(l + a)\p 

For LighthilPs gas a^ is not required; hence, 


- a 3 

v = a“ 
a l 


and 


a = 


a 5 


and in a manner similar to that given previously 


a r 


u = 


a l 


_ 2(7 + q)(l + q) f„ 
V a l/ (4 + a)' 


v2' 


~ ~ T t 


1/2 


7 + a 
4 + a 


(32) 

(33) 

(34) 


(35) 


(36) 


(for Lighthill gas). Then p, p, and T may be determined from equations (32), (33), 
and (34). 


Body Points 

As previously mentioned, the finite-difference scheme of equation (17) can be used 
for the body point n+1,0 only if imaginary flow properties are known at the point n,-l 
inside the body. Since the computational procedure cannot determine the properties at 
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k = -1, these properties must be determined by an extrapolation scheme or reflection 
principle. Several extrapolation and reflection techniques (such as those in ref. 13) were 
tested in the computation of the frozen flow past a wedge (for which the exact solution is 
known) by the present method. However, the results showed that the properties computed 
at the body points by all these techniques were less accurate than those computed by the 
method given herein. 



It was found that more accurate results were obtained when the properties at mesh 

points n-1,0 and n,l were used to determine the properties at the body point n+1,0 

as illustrated in figure 4. The finite-difference scheme replaces i PM by the sym- 

\ ^ j n 0 

metrical difference quotient ’ 

(Mn+l.O ~ ( am )n-l,0 
2 Ax 


replaces 



by the modified forward difference quotient 



n+1,0 + 

Ay 



and replaces (d m ) n Q by the average value 



( dm )n-l,0 
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With these replacements, the system of partial differential equations 


become 



9b. 


m 


3y 


+ d 


m 


n,0 


= 0 


(m = 1, 2, 3, 4, 5) 


( a m + = ( a m 


Axd 


m )n-l,( 


2 Ax 
Ay 




(37) 


for m = 1, 2, 3, 4, 5. 

At first glance equation (37) appears to be insufficient because (d m ) n+ i o an< * 

(bm) n+ i o a PP ear > and the properties at n+1,0 are unknown. However, because of the 

coordinate system used, there are two simplifications which allow a solution to be 
obtained. First, the m = 3 vector is not required because a 3 = puvri and the boundary 
condition on the body requires that v n+1 q = 0 for all values of n; hence, ( a ^=0. 
Second, since ’ 


bj = pvflx x 

b 2 = b x u 

b 4 = b lE v 


b 5 = bjQ! 


( b m) n+ i o = ( b m) n i o = 0 f° r m = 1, 2, 4, 5. With these simplifications, equation (37) 
becomes 


(a m + Axd m ) n+10 - (a m - Axd m ) n _^ 0 - (38) 

for m = 1, 2, 4, 5. 

For mesh points on the body, equation (38) gives (a m + Axd m ) n+ j q in terms of 

the known properties at n,l and n-1,0 for m = 1, 2, 4, 5. Also, it has previously 
been shown that ( a 3) n+ j q = 0 and (dl) n+ j g = 0. Explicit relations for the individual 

fluid properties can be obtained for a perfect (frozen) gas only, since d^ and dg are 
zero for these gases. 
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At the body surface, v = 0, y = 0, A x = 1, and d2 = -jp sin 9; therefore, 


Substitute equation (28) for 


a 2 + Axdg _ jj_/^ _ j Ax sin g\ + - 

a x pu\ " r / 


E into this equation to obtain 
2 

WjU + O^U + OJg = 0 


(3 


(4 


where 


u,-l- - L£Sjlnj\ 

1 7 + 3q!\ f J 

a2 + AXd2 

co 9 = l 

2 a l > 


co 3 = 


Ht - (1 - a)E v ip - aO D il/ 


2(1 + a) L _ j Ax sin 9 
7 + 3o! I r 


(4 


The solution of equation (40) for supersonic flow is 


_ ~ (t, 2 + ( a> 2 2 ~ 4^1^03] 


1/2 


2 001 


(4 


Equation (42) gives an explicit relation for u on the surface for a perfect (frozen) gas 
since 


a = Constant 


and 


E y = Constant 


for a frozen gas. The other surface variables may be determined now from 

a l 


P = 


ur^ 


- _ ( a 2 + ^2) - u(ai) 


r ^ - j Ax sin 9 



and 


T = E 

p(l + a)\p 


( 45 ) 


For the nonequilibrium gas, a and E v are not known and the quantities 

(a m + Axd m ) n+ j q for m = 1, 2, 4, and 5 cannot determine them explicitly. Therefore, 

an iterative technique must be used; one such technique is described in appendix B. 

To illustrate the accuracy of the present method compared with the reflection prin- 
ciple, the frozen flow over a 40.02° wedge was computed for y - 1.4 and M m = 6. The 
reflection principle sets 


u n,-l 


= u 


n,l 


P n,-1 = V 
p n,-l = p n,l 
v n,-l = -v n,l 


and uses equation (17) on the body 
surface. Figure 5 shows the density 
along the surface of the wedge. The 
superiority of the present method over 
the reflection principle is clearly evi- 
dent here. Actually, the reflection 
principle should be more accurate for 
a straight body than a curved body 
because it gives incorrect gradients 
normal to the surface of a curved body. 
(For instance, the reflection principle 

yields = 0 on the surface, but the 
8y 9n ’ 

centrifugal effects make + 0 on 

oy 

a curved surface.) 

For the Lighthill gas, the equa- 
tion analogous to equation (42) is 


5.2r 


5- 1 1— 


o 

A 


■ Exact solution 


o Present method 
a Reflection principle 


— a u ° a ^ 00 ooouoo 

P 5.0|~ A o o ° ° ° °°° 


4.9) 


4.81 


a a Aa a a AA a 


J L 


10 20 30 40 50 

Surface mesh point, n 


Figure 5.- Density along wedge surface for frozen flow. M 00 = 6; 
y = i.4; e = 40.02°. 


^2 ^^2 


u - 


+ Axd 


vH .\2 


- 4o) 4 co 5 


1/2 


2ci>/ 


(46) 
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where 


and 




_ i _ 1 + a L _ j Ax sin fl\ 

2(4 + a) \ r / 


- (^t ~ + L _ j Ax sin 0\ 

^5 4 + 0 ! \ r / 


( 47 ) 


The quantities p, p, and T may then be obtained from equations (43), (44), and (45), 
respectively. 

As shown in reference 22, the proper stability criterion for the body points, in 
supersonic flow, is that the n+1,0 mesh point should not fall beyond the intersection of 
the right-running frozen Mach line, from n,l, with the body surface. (See fig. 4.) Equa- 
tion (38) is consistent with equation (17) in that it is accurate to 0(Ax2). 


RESULTS 


The numerical technique developed herein is applied to frozen flow over a cone — 
parabolic-arc— cylinder, vibrational nonequilibrium flow over a wedge, and flow of 
Lighthill’s gas over a wedge. These results are compared with the methods of charac- 
teristics and integral relations. Additional results for coupled nonequilibrium vibration 
and dissociation are given in reference 16. All the numerical computations were made on 
an electronic data processing system. 

As has been mentioned, the present method is applicable only to flow regions where 
the local frozen Mach number is greater than unity. For the case of a pointed body, with 
an attached bow shock wave, the entire flow field is generally supersonic. The flow in 
the vicinity of the nose may be approximated by using the frozen flow properties or by 
using the tip gradients derived for two-dimensional bodies in references 3, 5, and 23. 

This approximation gives the (initial) flow properties required for the present method to 
compute the entire flow field. 

When the flow over a blunted body is considered, the subsonic-transonic region in 
the vicinity of the nose must first be determined by some other method, such as that of 
reference 4 or 6. Such results can then be used to establish the initial data line for the 
present method, the solution continuing downstream in the supersonic flow region over 
the body. 
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Frozen Flow Over a Cone — Parabolic-Arc — Cylinder 

This example illustrates the accuracy of the present method for supersonic flow 
(Moo = 2, y = 1.4) over an axially symmetric body with a discontinuous curvature, but 
having a continuous slope. Figure 6 shows the body shape. The cone— parabolic-arc 
junction lies at % = 0.05, and the parabolic-arc — cylinder junction lies at \ = 0.5. Note 
that the tip of the cone is not at | = 0, but at £ = -0.002777. 

The initial data line 
was chosen to be the line 
normal to the surface at the 
cone— parabolic-arc junc- 
tion. Since the flow is still 
conical along this line, the 
Taylor-Maccoll cone solu- 
tion (ref. 24) could be used 
for the flow properties. 

Three mesh points were 
chosen inside the shock 
layer: k = 0 on the surface, 
k = 2 midway between the 
surface and the shock wave, 
and k = 4 just aft of the 
shock wave. The other 
points on the initial data line 
(k = 6, 8, 10, . . .) were in 
the undisturbed free stream. 

With Ay = 0.0045946. 

Ax = 0.00459 for \ i 0.5, 

and Ax = 0.0055 for f> 0.5, the flow field was computed as far as £ = 1.05. Since Ay 
remains constant and the distance between the body and the shock wave increases with 
the number of mesh points inside the shock layer becomes much greater than three as the 
computations proceed downstream. The time required by the computing machine for this 
flow field was 2.2 minutes. 

The pressure coefficient along the body surface is plotted in figure 6 along with that 
given by the (standard) method of characteristics. This figure shows good agreement 
between the two methods. 



Figure 6.- Variation of surface pressure coefficient with axial distance for frozen flow past a 


cone— parabolic-arc— cylinder. = 2; y = 1.4; C* 


p V 2 

r oo oo 
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Nonequilibrium Flow Over Wedges 

It is well-known that for the supersonic, inviscid flow of a frozen gas past a wedge, 
the attached shock wave is straight and the fluid properties are constant between the body 
surface and the wave. The same qualitative results also hold for a gas in equilibrium; 
however, in this case, the shock wave is closer to the body and the constant fluid proper- 
ties are different from those of frozen flow. In particular, the equilibrium pressure and 
temperature are less than the corresponding properties for frozen flow. 

When the conditions are such that nonequilibrium effects must be considered, the 
flow field is complicated. At the tip of the wedge the flow is frozen since it has just 
passed through the infinitesimally thin shock wave and has not had time to relax. The 
translational and rotational energy modes reach equilibrium almost immediately; however 
the vibrational and dissociational energies require a finite time to reach equilibrium. 

For an undissociated free stream, the vibrational energy and the degree of dissociation 
increase along the surface of the wedge. These conditions in turn decrease the tempera- 
ture and pressure, but increase the velocity and density. The shock wave is inclined at 
the angle for frozen flow at the tip, but bends downward and approaches the equilibrium 
angle far downstream. As the flow along the surface achieves an equilibrium condition 
downstream, the pressure is the same as that for the equilibrium wedge solution; how- 
ever, the same is not true for the other flow variables (ref. 1). Downstream of the nose, 
three regions are observed in the profiles of properties between the body and the shock 
wave. Near the surface there exists an entropy layer since the streamlines passed 
through the steeper shock wave near the tip of the wedge. In the vicinity of the shock 
wave, the flow properties are close to being frozen since they have just passed through 
the shock wave and have not had time to relax. Between these two regions is a region 
which is essentially in equilibrium. This portion of the flow field has passed through the 
shock wave further upstream and thus has had time to relax to the equilibrium conditions. 

The wedge was chosen, both for simplicity and for comparison with other methods, 
in order to investigate the use of the present method for nonequilibrium flows. For the 
initial data line, a line normal to the surface and at a finite distance from the tip has been 
used for each case. The flow properties along this line have been determined (approxi- 
mately) by using the exact wedge-tip gradients given in appendix C. 

It was found that an initial data line with only two points between the wedge and the 
shock wave produced good results. In each case these two points were chosen as k = 0 
at y/6 = 0, and k = 2 at y/6 = 2/3. The remainder of the initial data line 
(k = 4, 6, 8, . . .) values were free-stream quantities. The distance of the initial data line 
from the tip x Q determines the spacing Ay, whereas Ax is limited by the stability 
criterion. Care must be exercised to insure that Ax remains small compared with the 
relaxation lengths. 
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Since the computations extend through the shock wave, the free stream, as well as 
the shock layer, must be considered in determining the stability criterion. In all the 
examples considered, the free- stream stability criterion was more stringent than that of 
the shock layer. This difference is due to the fact that the free -stream velocity is 
inclined at the wedge half-angle 0 to the coordinate system. 

As noted by Lax (ref. 9) and demonstrated in appendix A, the shock waves are nar- 
rowest when Ax is chosen as large as possible but still satisfies the stability criterion. 
Therefore, in the examples given, the largest value of Ax consistent with the stability 
criterion was used. It was found that this value was always small compared with the 
relaxation lengths of vibration and dissociation. 

Vibrational nonequilibrium .- Sedney et al. (ref. 1) and South (ref. 5) have solved the 
vibrational nonequilibrium flow field over a wedge by the methods of characteristics and 
integral relations, respectively. The following case (designated case I) was chosen to 
compare with the method of characteristics: 

Case I: The gas is N 2 and 


T m = 300° K 

Moo = 6 

n.12 

Too 


e = 40.02° 
Woo = 0 


All lengths were nondimensionalized by 


Therefore, the proportionality constant in the relation for r does not affect the nondi- 
mensional variables. 

In order to determine the effect of the mesh spacing, two separate computations were 
performed with initial data lines at x Q = 0.035 and at x 0 = 0.1* Although Ax and Ay 
were smaller for the case with x Q = 0.035 than those for x 0 = 0.1> the ratio Ax/Ay 
was the same for both cases. The flow field was computed downstream to x = 4 in each 
case. Figures 7, 8, 9, and 10 show the variation of p, T, E v , and u, respectively, 
along the wedge surface. Agreement with the method of characteristics is very good, 
except near the starting point. Here the pressure and velocity, in particular, deviate from 
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Figure 7.- Pressure along wedge surface for case I. 

the correct solution; however, these quantities quickly recover as the calculations are 
made downstream. The amount of deviation increases with increasing mesh size. On the 
other hand, figures 8 and 9 show very little deviation in the values of T and E y from 
the characteristics solution. 

Figures 11 and 12(a) give the p and T profiles normal to the surface at x = 4. 
Of particular interest are the results in the vicinity of the shock wave where the values 
of p and T, in the shock layer, drop to free-stream values over a few mesh points; 
thus the shock wave is smeared over these points. Also, these figures show that a more 
definitive shock is produced by decreasing the mesh spacing (while Ax/ Ay is held 
constant). (It is interesting to note that the computational time was 8.8 minutes for 
x Q = 0.035 and only 2.0 minutes for x Q = 0.1.) 
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Figure 12(b) shows a portion of the temperature profile at an expanded scale. This 
figure illustrates the entropy layer near the surface, a region of near -equilibrium, and the 
nonequilibrium region near the shock wave proper. Also, there is a significant difference 
between the actual shock location and that for frozen flow. This difference is due to the 
fact that the shock wave has the frozen slope at the tip, but bends downward and approaches 
the slope of an equilibrium flow as x increases. (See ref. 1.) 

Lighthill's ideal dissociating gas .- Capiaux and Washington (ref. 3) have used the 
method of characteristics to compute the nonequilibrium flow field over a wedge for 
Lighthill's ideal dissociating gas. Their technique, however, differs from that used by 
Sedney in that they used the stream function and distance perpendicular to the wedge axis 
as independent variables. Newman (ref. 7) also solved this problem by using a modified 
method of integral relations. The following case (designated case n) was computed by the 
present method and corresponds to case I in reference 3: 
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Figure 10.- Velocity along wedge surface for case I. 
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Figure 13.- Variation of surface pressure with axiai distance for case 1 1. 
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Figure 14.- Variation of surface temperature with axial distance for case 1 1. 



All lengths were nondimensionalized by the same length scale as that used in reference 3, 
that is, 


L = 


~oo 


2.5 


In this manner the numerical value of C* does not affect the nondimensional properties. 
Note, however, that reference 7 used a different length scale. 

In order to obtain good accuracy near the tip of the wedge, the problem was solved 
in two parts. First, with an initial data line at x 0 = 10"^, the solution was carried out to 
x = 4 x 10“3. By using the results already computed at x = 1.18 x 10"^, a new initial data 
line (of two points in the shock layer) was used to carry the solution to x = 0.111. The 
second solution partly overlapped the first; therefore, the fluctuations which occur near the 
starting point would not be present in the region 4 x 10"^ = x ^ 0.111. These results are 
plotted in figures 13, 14, and 15 and show the variation of p, T, and a, respectively, 
along the surface. It should be noted that reference 3 measures x along the axis of 
symmetry, whereas the present x is along the surface. Therefore, x cos 9 in this 
report corresponds to x in reference 3. 

Figure 13 shows that the deviation of p, near the starting point, is more pronounced 
than in the previous example for vibrational nonequilibrium. However, it should be noted 
that x cos 9 is plotted on a logarithmic scale; hence, this deviation extends over a very 
small distance on a linear scale. For x cos 9 > 2 x 10 " 3 , the present method and that of 
references 3 and 7 all disagree on the surface pressure. It is not presently known why 
these results differ, or which is the more correct solution. 

Figure 14 shows that the present method and that of reference 7 give the same 
results for T along the surface. The method of reference 3, however, differs from 
these two methods for x cos 9 > 3 x 10“^. All three methods give essentially the same 
results for a., as shown in figure 15. (The total computational time for the present 
method was 16.2 minutes.) 


DISCUSSION 

The comparisons given in figures 5 to 15 show that the present method correctly 
calculates the flow field over bodies except near the initial data line. At that point, 
unrealistic deviations appear in the flow properties on the body but are quickly damped 
out. It was found that these deviations still occurred even when as many as 20 points were 
used on the initial data line between the body and the shock wave. The only difference was 
that the deviations started at a position downstream where the first mesh point outside the 
shock wave on the initial data line influences a body point. It is believed that these 
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deviations are caused by using free-stream properties at one mesh point, and shock prop- 
erties at the mesh point adjacent to it on the initial data line. This procedure causes a 
weak disturbance to propagate down to the body, and the boundary conditions cause the sur- 
face properties to deviate from the correct solution. When other types of coordinate sys- 
tems and schemes for computing the body points were investigated, the deviations near the 
initial data line were greater than those of the present method. Also, it was mentioned 
previously that the y-momentum equation was not used at the surface because it is not 
needed. (The boundary condition specifies that v = 0 on the surface.) Upon investi- 
gating the computed flow fields, it was found that the y-momentum equation was not satis- 
fied at the surface in the regions where these deviations occurred. However, downstream 
where the deviations were damped out, it was found that the y-momentum equation was 
satisfied at the surface. The present method could be altered to satisfy simultaneously 
all the equations and boundary conditions on the surface if Ax was made a variable for 
each step in the computational procedure. However, this procedure increases the com- 
plexity of the method manyfold and may require computational times comparable with the 
method of characteristics. 

Several investigations were performed to determine the effect of using different 
values of the ratio Ax/ Ay. When this ratio was greater than that given by the stability 
criterion, the computed flow properties became unstable after only a few mesh spaces 
downstream of the initial data line. However, for values of Ax/Ay consistent with the 
stability criterion, more accurate results were always obtained the closer the mesh 
system followed the characteristics. This result is illustrated in appendix A. 

The staggered mesh system that is used for the finite-difference scheme has the 
decided advantage that it requires only half as many mesh points as an unstaggered sys- 
tem and yet it gives the same accuracy. This system reduces both the computing machine 
storage spaces and the computational time by a factor of 2. The present method also has 
the advantage of using a mesh system that is fixed, whereas the method of characteristics 
must compute its mesh points. Therefore, the present method requires less computa- 
tional time than the method of characteristics. 

The results obtained show that the present method gives good accuracy even for a 
large mesh scale. For the flow past a wedge, only two points were used inside the shock 
layer for the initial data line. The downstream profiles of the flow properties across the 
shock layer show that the present method gives better accuracy than the method of integral 
relations. If greater accuracy is desired in the present method, an extrapolation proce- 
dure such as described by Roberts (ref. 10) can be used. 

In addition to the frozen flow past axially symmetric bodies, the present method has 
also been applied to nonequilibrium flows over cones. However, since there are no exact 
expressions for the cone tip gradients like those for wedges, the initial data line was 
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composed of the frozen cone solution. As a result, the solutions were not very accurate 
for the mesh spacings used to date. Reference 2 also encountered considerable difficulty 
in the computation of nonequilibrium flows over cones by the method of characteristics. 

In reference 2 it was found that a much smaller mesh spacing was needed for cones than 
for wedges. It is believed that accurate, nonequilibrium cone solutions can be obtained 
with the present method if a very small mesh spacing was used, or if the cone tip gradi- 
ents were known. 



Figure 15.- Variation ol surface atom mass fraction with axial distance for case 1 1. 

CONCLUSIONS 

The finite-difference method developed by P. D. Lax for unsteady flows is applied 
to the steady-state equations for the supersonic region of inviscid flows past two- 
dimensional and axisymmetric bodies. A diatomic gas subject to nonequilibrium vibra- 
tions and dissociation is considered. As a consequence of Lax's finite -difference scheme, 
an artificial viscosity is implicitly introduced. This scheme allows the computations to 
proceed downstream of an initial data line as if no shock wave were present at all. The 
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shock wave appears in the solution, however, smeared over several mesh spaces while it 
accurately gives the proper jump conditions across the shock. As a result of the applica- 
tion of this method, the following conclusions can be drawn: 

1. The finite-difference method developed herein is applicable to the solution of the 
supersonic -hypersonic flow fields past two-dimensional and axisymmetric bodies. 

2. The numerical results show that the present method compares favorably with the 
results obtained by the methods of characteristics and integral relations for both perfect- 
and real- gas flow fields. 

3. The new difference scheme developed for the body points is more accurate than 
extrapolation schemes and reflection principles previously used. The success of this 
difference scheme is dependent on a body-oriented coordinate system. 

4. The staggered mesh system has the advantage of reducing the computations 
required for an unstaggered mesh system without sacrificing accuracy in the results. 

5. The closer the mesh system approaches the characteristics network, the more 
accurate the results. 

6. The present method is sufficiently accurate to yield better profiles of properties 
across the shock layer than the method of integral relations, yet it is simple enough to 
require less time for machine computations than the method of characteristics. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., September 21, 1965. 
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APPENDIX A 


FROZEN FLOW THROUGH AN OBLIQUE SHOCK WAVE 

This example illustrates the application of the present method for the computation 
of a steady flow field without a boundary and is similar to the computations by Lax 
(ref. 9) for unsteady flow. This case is for an oblique shock wave (ftp = 30°) positioned 
in a free stream of = 3 and y = 1.4. The x-axis is parallel to V^, and the y-axis 
is perpendicular to V M . Since there is no boundary, the initial data line (n = 0) consists 
of the free-stream properties at k = 1, 3, 5, 7, . . . and the frozen shock properties at 
k = -1, -3, -5, -7, . . . and places the shock wave at k = 0 on the initial data line. The 



Figure 16,- Pressure profile through an oblique shock wave at n - 49 for frozen flow. 
M oo = 3: ° F = 30 ° : r = 1A 
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APPENDIX A 


present method has been used to compute the downstream flow field for three values of 
Ax/Ay. From the stability criterion, 


=£ = cot tan' 1 a +p 

AY /max Wf 


= 1.2910 


Figure 16 shows the pressure profiles computed at n = 49. This figure illustrates that 
the profile through the shock wave is qualitatively similar to those obtained experimen- 
tally in reference 19. It is also seen that the computed shock wave is narrowest when 
= (x£) > a conclusion noted by Lax in reference 9. The reason for this phenomenon 

is that the "artificial" coefficient of viscosity I i n eq. (19)1 decreases as Ax/Ay 
increases (Ay being held constant) . 
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APPENDIX B 


BODY POINTS FOR NONEQUILIBRIUM FLOW 


The numerical scheme formulated for the body points computed a m + Axd m for 
m = 1, 2, 4, and 5 at the new mesh point n+1,0 from the known properties at n,l and 
n-1,0. For nonequilibrium flow a and E v are not known at n+1,0 and, hence, an 
iterative technique must be employed to determine the individual fluid properties from 
(a m + Axd m ^ n+ j q. The steps in this method are: 


(1) Assume 


( d4 )n+l,0 " Mn-1,0^ 


r n+l,0\ j 

0 / 


and 


( 2 ) 


( d 5)„ + l,o = ( d 5)„- 1 ,o(^) 

Compute (E v ) n+10 and o n+10 from 


( E v) 


n+1,0 


( a 4 + ^In+l.O (t»4) ntl ,o & 


( a l) 


n+1,0 


( a i) 


n+1,0 


and 


_ ( a 5 + ( d 5)n+i 0 

^ +1 ’° R 


Ax 


' n+1,0 


( a l) 


n+1,0 


(3) Compute u n+10 from equation (42), p n+1 Q from equation (43), p n+1 0 from 
equation (44), and T n+ j g from equation (45). 

(4) By using the properties computed in steps (2) and (3), compute new values of 
( d 4) n +l,0 and ( d 5) n +l,0 from equation (16). 

(5) By using the new values of (d4) n+ i o and ( d s) n +i g> repeat steps (2) to (4) until 
the temperature found in step (3) differs from its previous value by less than 10" ^ percent. 
Convergence usually requires less than five iterations. 
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APPENDIX C 


WEDGE -TIP GRADIENTS 


Exact expressions for the wedge-tip gradients have been derived by Sedney et al. 
(ref. 1) and South (ref. 5) for vibrational nonequilibrium, and by Capiaux and Washington 
(ref. 3) for Lighthill's ideal dissociating gas. South's derivation gives the results in a 
form suitable for obtaining the initial data line in the present method. Therefore, the 
results given in this appendix were obtained in the same manner as outlined in appen- 
dixes A and C of reference 5 but were extended to include dissociation. 


For frozen shock-wave relations, 

= l±3o k 
00 5 + 

Define F and G as 


then 


F = M oo 2 sin 2 n F 
c _ 2(F - 1) 


cot 9 = tan 0 F I — - 1 


u F = (1 - G)cos 9 + G cot sin 9 


(Cl) 

(C2) 

(C3) 

(C4) 

(C5) 


and 



G)sin 6 + G cot Ott cos 9 
£ 


- 2y»F - (r» - 1) 

F (y + l)y M 2 

\< OO / » OO OO 


PjT = 


fea + i ) f 

(y„ - ') F + 2 


(C6) 

(C7) 

(C8) 



- froo - ^ + 2 ] 

(^OO + 


(C9) 
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When derivatives of equations (C5), (C6), and (C7) are taken with respect to f2 F , 

fdu\ = ~ 4 cos % sin ( fl F - d ) G sin 6 


(^oo + X ) 


sin^fi-. 


dv\ _ 4 cos ftp cos - o) G cos q 
k dfi/F (y__ + 1) sin 2 ftp 


(^oo + X ) 

d p^ 4 sin fip cos Op 
, d ^/F (t~ + 1) 


(CIO) 


(C11) 


(Cl 2) 


With 6 defined as the y-position of the shock wave and 6 = the wedge-tip gradients 

L 

are 

- 2 (^~ "• ~ e O ~ 4T j 

dn . T f (7 + 3 a ,) " T r (7 + 3aJ (1 + «„) 


dx 




(Cl 3) 


duk 

dx 
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APPENDIX C 


By using these tip gradients and neglecting second and higher order terms, the flow prop- 
erties across the shock layer at the distance x Q (n = 0) from the tip may be written as 


Q o,k 


u. 



(C19) 


and 



(C20) 


P 0 ,k = PF + 




(C 21) 


(C22) 


“o,k 


= “oo + 


do'j 

dx 


(C23) 


By using these properties, T ■ may be obtained from the energy equation and then 
P 0 from the equation of state. Since = tan(S2 - 9 ), the shock position 6 at x Q 
is (third and higher order terms being neglected) 


5 = tan^p - o)x Q + sec^fip - 0^ 


df2\ x o 2 


dx/ 2 


(C24) 


When Lighthill's gas is considered, the only changes required are that equations (Cl) and 
(C13) be replaced by 
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oo 


4 + Ooo 

3 


(C25) 


dft 

dx 



and equations (Cl 7) and (C22) are not needed. 


(C26) 
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